{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_occupancy(res, phase, p_map, t_map):\n",
    "    r'''\n",
    "    Custom function to take the results of the percolation algorithm and update phase occupancy\n",
    "    '''\n",
    "    phase['pore.occupancy'] = 1 - res['pore.occupancy'][p_map]\n",
    "    phase['throat.occupancy'] = 1 - res['throat.occupancy'][t_map]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _calc_eff_prop(self):\n",
    "    r\"\"\"\n",
    "    Returns the main parameters for calculating the effective property\n",
    "    in a linear transport equation.  It also checks for the proper\n",
    "    boundary conditions, inlets and outlets.\n",
    "    \"\"\"\n",
    "    import numpy as np\n",
    "    if self.settings['quantity'] not in self.keys():\n",
    "        raise Exception('The algorithm has not been run yet. Cannot ' +\n",
    "                        'calculate effective property.')\n",
    "\n",
    "    # Determine boundary conditions by analyzing algorithm object\n",
    "    inlets, outlets = self._get_inlets_and_outlets()\n",
    "    Ps = np.isfinite(self['pore.bc_value'])\n",
    "    BCs = np.unique(self['pore.bc_value'][Ps])\n",
    "    Dx = np.abs(np.diff(BCs))\n",
    "    net = self.project.network\n",
    "    # Fetch area and length of domain\n",
    "    [amax, bmax, cmax] = np.max(net['pore.coords'], axis=0)\n",
    "    [amin, bmin, cmin] = np.min(net['pore.coords'], axis=0)\n",
    "    lx = amax-amin\n",
    "    ly = bmax-bmin\n",
    "    lz = cmax-cmin\n",
    "    A = lx*ly\n",
    "    L = lz\n",
    "    flow = self.rate(pores=inlets)\n",
    "    D = np.sum(flow)*L/A/Dx\n",
    "    return D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bulk_diffusion_wu(physics,\n",
    "                      network,\n",
    "                      phase,\n",
    "                      geometry,\n",
    "                      propname,\n",
    "                      diffusivity = 'pore.diffusivity',\n",
    "                      molar_density = 'pore.molar_density',\n",
    "                      throat_diameter = 'throat.diameter',\n",
    "                      throat_length = 'throat.length',\n",
    "                      pore_diameter = 'pore.diameter',\n",
    "                      **params):\n",
    "        r\"\"\"\n",
    "        Calculate the diffusive conductance of throats in network (instead of a\n",
    "        conduit) based on the areas\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        network : OpenPNM Network Object\n",
    "\n",
    "        phase : OpenPNM Phase Object\n",
    "        The phase of interest\n",
    "\n",
    "        Notes\n",
    "        -----\n",
    "        This function requires that all the necessary phase properties already be\n",
    "        calculated.\n",
    "\n",
    "        \"\"\"\n",
    "        #ct = phase.get_data(prop='molar_density',throats='all',mode='interpolate')\n",
    "        #Interpolate pore values to throats\n",
    "        DABt = phase.interpolate_data(propname='pore.diffusivity')\n",
    "        #Find g for full throat\n",
    "        tdia = network[throat_diameter]\n",
    "        tlen = network[throat_length]\n",
    "        gt = (np.pi*DABt*tdia**2)/(tlen*4)\n",
    "        g = gt[geometry.throats()]\n",
    "        phase[propname]=g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulation(n=8, npts=51):\n",
    "    r'''\n",
    "    Run a percolation and diffusion simulation\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    n : int\n",
    "        number of pores along each direction of the network used for diffusion.\n",
    "        Percolation network has length 2n along the z-axis\n",
    "\n",
    "    npts : int\n",
    "        number of percolation points\n",
    "    '''\n",
    "    import openpnm as op\n",
    "    import numpy as np\n",
    "    print('-' * 80)\n",
    "    print('Running Sim with Net Size =', n, 'num points', npts)\n",
    "    print('-' * 80)\n",
    "    Lc = 25e-6\n",
    "    pn1 = op.network.Cubic(shape=[n,n,2*n], spacing=Lc)\n",
    "    #code to run if boundaries was set to false\n",
    "    Ps = pn1.pores()\n",
    "    Ts = pn1.throats()\n",
    "    geo1 = op.geometry.GenericGeometry(network=pn1, pores=Ps,throats=Ts)\n",
    "    \n",
    "    low = .5e-6\n",
    "    high = 9.5e-6\n",
    "    \n",
    "    geo1['pore.diameter'] = 24e-6\n",
    "    radii = low + np.random.rand(pn1.num_throats())*(high - low)\n",
    "    geo1['throat.diameter'] = radii*2\n",
    "    geo1.add_model(propname='pore.volume', model=gm.pore_volume.sphere)\n",
    "    geo1.add_model(propname='pore.area', model=gm.pore_cross_sectional_area.sphere)\n",
    "    geo1.add_model(propname='throat.length', model=gm.throat_length.classic)\n",
    "    geo1.add_model(propname='throat.volume', model=gm.throat_volume.cylinder)\n",
    "    geo1.add_model(propname='throat.area', model=gm.throat_cross_sectional_area.cylinder)\n",
    "    \n",
    "\n",
    "    water = op.phases.Water(network = pn1)\n",
    "    Ps = geo1.pores()\n",
    "    Ts = geo1.throats()\n",
    "    phys_water = op.physics.Standard(network=pn1, phase=water, geometry=geo1)\n",
    "    IP_1 = op.algorithms.InvasionPercolation(network=pn1)\n",
    "    IP_1.setup(phase=water)\n",
    "    inlets = pn1.pores('bottom')\n",
    "    IP_1.set_inlets(pores=inlets)\n",
    "    IP_1.run()\n",
    "\n",
    "    z_vals = np.unique(pn1['pore.coords'][:, 2])[int(n/2):int(3*n/2)]\n",
    "    pore_map = np.in1d(pn1['pore.coords'][:, 2], z_vals)\n",
    "    throat_map = pn1.find_neighbor_throats(pn1.pores()[pore_map], mode='intersection')\n",
    "    pn2 = op.network.Cubic(shape=[n,n,n], spacing=Lc)\n",
    "    prj2 = pn2.project\n",
    "    Ps = pn2.pores()\n",
    "    Ts = pn2.throats()\n",
    "    geo2 = op.geometry.GenericGeometry(network=pn2, pores=Ps, throats=Ts)\n",
    "    for prop in geo1.props():\n",
    "        if 'pore' in prop:\n",
    "            geo2[prop] = geo1[prop][pore_map]\n",
    "        else:\n",
    "            geo2[prop] = geo1[prop][throat_map]\n",
    "        \n",
    "    air = op.phases.Air(network = pn2)\n",
    "    phys_air = op.physics.Standard(network=pn2, phase=air, geometry=geo2)\n",
    "    del phys_air['throat.diffusive_conductance']\n",
    "    bulk_diffusion_wu(physics = phys_air, network=pn2, phase=air, geometry=geo2, propname='throat.diffusive_conductance')\n",
    "    phys_water.regenerate_models()\n",
    "    phys_air.regenerate_models()\n",
    "\n",
    "    x_values = []\n",
    "    y_values = []\n",
    "    \n",
    "    res = IP_1.results(0.0)\n",
    "    update_occupancy(res, air, pore_map, throat_map)\n",
    "    phys_air.add_model(model=op.models.physics.multiphase.conduit_conductance,\n",
    "                       propname='throat.conduit_diffusive_conductance',\n",
    "                       throat_conductance='throat.diffusive_conductance',\n",
    "                       mode='strict')\n",
    "    bulk_diffusivity = air['pore.diffusivity']\n",
    "    bottom_boundary = pn2['pore.bottom']\n",
    "    top_boundary = pn2['pore.top']\n",
    "\n",
    "    for x in range(npts):\n",
    "        IPsat = float(x)/float(npts)\n",
    "        res = IP_1.results(IPsat)\n",
    "        update_occupancy(res, air, pore_map, throat_map)\n",
    "        phys_air.regenerate_models(propnames=['throat.conduit_diffusive_conductance'])\n",
    "        mask = air['throat.conduit_diffusive_conductance'] == air['throat.diffusive_conductance']\n",
    "        am = pn2.create_adjacency_matrix(weights=mask, fmt='coo')\n",
    "        if tt.ispercolating(am,\n",
    "                            bottom_boundary,\n",
    "                            top_boundary,\n",
    "                            'bond'):\n",
    "            print('Step', x, 'is percolating')\n",
    "            Fickian_alg = op.algorithms.FickianDiffusion(network=pn2, phase=air)\n",
    "            Fickian_alg.set_value_BC(values=0.6, pores=top_boundary)\n",
    "            Fickian_alg.set_value_BC(values=0.2, pores=bottom_boundary)\n",
    "            Fickian_alg.setup(conductance='throat.conduit_diffusive_conductance')\n",
    "            Fickian_alg.run()\n",
    "            effective_diffusivity = _calc_eff_prop(Fickian_alg)/np.mean(air['pore.molar_density'])\n",
    "            prj2.purge_object(Fickian_alg)\n",
    "        else:\n",
    "            print('Step', x, 'is NOT percolating')\n",
    "            effective_diffusivity = 1e-12\n",
    "        #calculation of saturation\n",
    "        p_vol = pn2['pore.volume']\n",
    "        saturation = 1 - np.sum(air['pore.occupancy']*p_vol)/np.sum(p_vol)\n",
    "        Deff = (effective_diffusivity/bulk_diffusivity)[0]\n",
    "        print('Step', x, 'Saturation', \"%.3f\" % saturation, 'Diffusivity', \"%.3f\" % Deff)\n",
    "        x_values.append(saturation)\n",
    "        y_values.append(Deff)\n",
    "\n",
    "    return x_values, y_values"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
